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We present a numerical code designed to study astrophysical phenomena involving dynamical 
spacetimes containing black holes in the presence of relativistic hydrodynamic matter. We present 
evolutions of the collapse of a fluid star from the onset of collapse to the settling of the resulting 
black hole to a final stationary state. In order to evolve stably after the black hole forms, we excise 
a region inside the hole before a singularity is encountered. This excision region is introduced after 
the appearance of an apparent horizon, but while a significant amount of matter remains outside 
the hole. We test our code by evolving accurately a vacuum Schwarzschild black hole, a relativistic 
Bondi accretion flow onto a black hole, Oppenheimer-Snyder dust collapse, and the collapse of 
nonrotating and rotating stars. These systems are tracked reliably for hundreds of M following 
excision, where M is the mass of the black hole. We perform these tests both in axisymmetry and in 
full 3+1 dimensions. We then apply our code to study the effect of the stellar spin parameter J/M 2 
on the final outcome of gravitational collapse of rapidly rotating n = 1 polytropes. We find that a 
black hole forms only if J/M 2 < 1, in agreement with previous simulations. When J/M 2 > 1, the 
collapsing star forms a torus which fragments into nonaxisymmetric clumps, capable of generating 
appreciable "splash" gravitational radiation. 
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I. INTRODUCTION 

Since many of the most interesting phenomena in as- 
trophysics involve black holes, the modeling of black hole 
spacetimes is one the most important problems in numer- 
ical general relativity. It is also one of the most challeng- 
ing problems. Black hole evolutions present all the usual 
difficulties of numerical relativity, such as the need to find 
a stable form of the field evolution equations and the need 
to find a practical coordinate system. In addition, han- 
dling the singular region is very subtle for a numerical 
code; the black hole singularity must be avoided to allow 
the exterior evolution to continue far into the future. 

One of the most promising methods to date of dealing 
with black hole singularities is black hole excision. This 
method, first suggested by Unruh 0, exploits the fact 
that the singularity resides inside an event horizon, a re- 
gion that is casually disconnected from the rest of the 
universe. Since no physical information propagates from 
inside the event horizon to outside, one should be able 
to evolve the exterior independent of the interior space- 
time. Inside the event horizon, causality entitles us to do 
anything which will produce a stable exterior evolution. 
In particular, one can excise a region inside the horizon 
containing the singularity and replace it with suitable 
boundary conditions at its outer surface. 

Although it is guaranteed that no physical signal can 
propagate from inside the horizon to outside, unphysical 
signals often can propagate in evolution codes. Gauge 
modes can move acausally for many gauge conditions. 
Although they carry no physical content, such modes 
may destabilize the code. Thus, the choice of gauge is 
crucial to obtaining good excision evolutions. In addi- 
tion, constraint-violating modes can, for some formula- 



tions of the field equations, propagate acausally, creating 
inaccuracies and instabilities. Thus, the choice of for- 
mulation is also crucial to obtaining good excision evolu- 
tions. 

The feasibility of black hole excision was demonstrated 
in spherically symmetric 1+1 dimensional evolutions of 
a single black hole in the presence of a self-gravitating 
scalar field 0- 0- 0- Excision was also implemented 
successfully to study the spherically symmetric collapse 
of collisionlcss matter to a black hole in Brans-Dicke the- 
ory [6j . Three-dimensional evolutions of black holes with 
excision were attempted by using the standard 3+1 ADM 
formulation, for a stationary 7] and for a boosted black 
hole [8| . Although the introduction of excision improved 
the behavior of these black hole simulations, long-term 
stability could not be achieved due to instabilities en- 
demic to the unmodified ADM formulation. 

Since then, new and more stable formulations of the 
3+1 Einstein field equations have been devised. Using 
excision in a modified version of the ADM equations com- 
monly referred to as BSSN HG3, 

several groups have 
evolved stationary black hole spacetimes (non-spinning 
and spinning) for arbitrarily long times [TH Il2j . Long- 
term stability has also been achieved usin g hy perbolic 
formulations of the field equations [l3L ItS, Il5| and us- 
ing characteristic evolutions [lg. Success has also been 
achieved in evolving distorted and moving black holes 
with excision, both with characteristic formalisms fl6| 
and with BSSN [llE[l<|. Excision has also been used 
to simulate the grazing collision of two black holes [2(j 
and to simulate binary black holes for approximately one 
orbital period [2l| . 

The last several years also have seen significant ad- 
vances in numerical, 3+1 relativistic hydrodynamics in 
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dynamical spacetimes (see, e.g. ) • The simu- 

lation of rapidly rotating, relativistic stars is now pos- 
sible, and fully relativistic evolution codes are being 
used to study the stability [2(| and gravitational col- 
lapse jU H3, EJ of such objects. Binary neutron stars 
can now be evolved accurately for multiple orbits 
and numerical simulations have been used to study the 
final merger of these binaries 31, 32]. 

The necessary tools are clearly being forged to enable 
numerical relativity to model a wide variety of strong- 
field gravitational phenomena. Many interesting systems 
in astrophysics involve the simultaneous presence of both 
black holes and hydrodynamic matter fields, and these 
systems will require a code which can handle both in 
order to model them reliably. 

One important scenario involving both hydrodynamic 
matter and black holes is core collapse in massive stars, 
an event of immense importance due to its association 
with supernovae, the formulation of neutron stars and 
black holes, and gamma ray bursts (GRBs). Some recent 
numerical simulations suggest |33| that a star must have 
mass less than about 20 Mq for core collapse to result in 
a conventional neutron star and supernova explosion. For 
progenitor masses between around 20M Q and 40M Q , the 
core collapses to a neutron star initially, but it eventually 
implodes to a black hole, as ejected material slowly falls 
back onto the remnant (see also [34[ ). For more mas- 
sive stars, the core collapses promptly to a black hole. 
Such a massive system is a promising candidate for a 
GRB 35| . There is growing evidence that long duration 
GRBs are associated with hypernovae that accompany 
the collapse of massive stellar cores. This evidence in- 
cludes the association of the low-energy GRB980425 with 
a supernova j3||, the presence of supernova- like features 
in the optical afterglow of several GRBs j33, and the 
existence of freshly synthesized elements in the ejecta of 
GRB 011211 [3^. Most recently, a hypernova was found 
to be temporally and spatially coincident with a normal 
cosmological burst source, GRB 030329 [3!j. Most mod- 
els of the central engine of GRBs involve a black hole 
surrounded by a rapidly accreting disk and a jet [4p| . 
Three dimensional fully relativistic simulations of both 
the black hole and the exterior matter will be needed 
to test the feasibility of various models for the produc- 
tion of GRBs from such "collapsars" , and it is likely that 
excision will be required to track the full evolution. 

The merger of binary neutron stars is a promising 
source of gravitational waves, as well as a prime can- 
didate for short duration GRBs 41J. Binary mergers of 
stars of high compaction collapse promptly to a black 
hole [U- The coalescence of low-compaction neutron 
stars probably leads to the formation of a hypermassive 
neutron star remnant [26l l32l 142^ followed by a delayed 
collapse [43L |44| . Either way, excision of the black hole 
singularities is necessary to follow binary mergers which 
form black holes with accretion disks, emit gravitational 
waves, and drive short-duration GRBs. 

The dynamics of accretion flows onto a black hole is 



another problem of great importance, since black holes 
are usually visible electromagnetically only through ac- 
cretion. When the mass of the accreting fluid is much 
less than that of the black hole, then the matter can 
be evolved on a fixed black hole spacetime. When the 
masses of the hole and the matter are comparable, then a 
fixed black hole spacetime becomes a bad approximation 
to the true metric, and the full system must be evolved 
self-consistently. This is particularly important for deter- 
mining if and when the disk may produce an instability, 
as in the runawa y in stability |45j, or in the one-armed 
spiral instability 1461 which can generate quasi-periodic 
gravitational waves |47| . 

Other coupled black hole-hydrodynamic matter sys- 
tems include a neutron-star black hole binary, which 
might be an important source of gravitational radiation 
for LIGO, and also the disruption or capture of a star 
by a massive black hole, which is expected to be a major 
source of waves for LISA [4£| • Supermassive black hole 
seed formation by the collapse of a massive or supermas- 
sive star is another important example |4*7Ll49| . 

Successful attempts at evolving matter and a black 
hole together in a dynamical spacetime using excision 
have been rare. Scheel el al. [fj simulated the collapse 
of a spherically symmetric configuration of collisionless 
matter in Brans-Dicke theory with excision. The space- 
time within the numerical domain was evolved until the 
appearance of an apparent horizon. At that time, an 
excision boundary was introduced and the evolution of 
the exterior spacetime was continued. An attempt to 
evolve a dynamical black hole spacetime with hydrody- 
namic matter was undertaken by Brandt el al [50j in ax- 
isymmctry. In that paper, a black hole is evolved with an 
accretion flow, using the ADM formalism and an isom- 
etry inner boundary condition at the apparent horizon. 
They were able to evolve several systems for up to about 
lOOAf . Little progress has been made since then, pre- 
sumably because the computational tools for performing 
excision and relativistic hydrodynamics in 3+1 had to be 
perfected independently. We have only now reached the 
stage where these tools can be put together successfully. 

In this paper, we perform the first simulations which 
utilize excision to evolve relativistic hydrodynamic mat- 
ter in 3+1 dynamical spacetimes containing black holes. 
In particular, we present evolutions of the gravitational 
collapse of stars from the beginning of collapse, through 
black hole formation, to quiescent final states. We per- 
form these evolutions in two stages. From the begin- 
ning of collapse until the appearance of an apparent hori- 
zon, we evolve using our new, relativistic hydrodynamics 
(BSSN) code without excision (i.e. our "pre-excision" 
code). After an apparent horizon appears, we continue 
the evolution with a region inside the horizon excised. At 
the moment we introduce excision, a significant amount 
of matter is still outside the excision zone and the black 
hole is significantly distorted and in a nonstationary 
state. We follow its evolution to a final stationary state. 
In Section[n]we describe our evolution scheme, including 
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our gauge and boundary conditions. In Section II I II we 
describe our code diagnostics. In Section HVI we test our 
code on systems with known behavior, including vacuum 
black holes, relativistic Bondi accretion, Oppenhcimer- 
Snyder collapse, and the collapse of unstable polytropes, 
both non-rotating and rotating. We find that we are able 
to evolve many systems stably for hundreds of M . When 
we evolve systems with appreciable angular momentum, 
we can only conserve J for somewhat shorter durations, 
but this duration can be extended by increasing reso- 
lution. In Section [V] we apply our code to study the 
late-time outcome of pressure-depletion-induced gravita- 
tional collapse of rapidly-rotating polytropes with poly- 
tropic index n = 1. We find that stars with J/M 2 < 1 
collapse to Kerr black holes with no surrounding disks. 
Stars with J/M 2 > 1 collapse to tori, which then frag- 
ment. This fragmentation process can produce copious 
amounts of gravitational radiation, originally referred to 
as "splash radiation" |51| . Finally, we summarize our 
results and discuss future improvements to our code in 
Section EH 

Throughout this paper, Latin and Greek indices denote 
spatial components (1-3) and spacetime components (0- 
3), respectively. We use geometrized units, so that G — 
c= 1. 



variables must satisfy the following constraint equations: 
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These constraints are solved only at the initial time, and 
are used henceforward as diagnostics. In an attempt to 
improve the stability and accuracy of the BSSN formu- 
lation, one can add multiples of the above constraints to 
the field equations. Many possible modifications of this 
kind have recently been suggested \56L I57L \58L l59l l60| . 
We found a slight improvement in ADM mass conserva- 
tion by adopting the following modifications: 



d t cf> = 
d t % = 
d t Ai 3 = 
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(8) 



II. SUMMARY OF METHOD 



Our basic code has been described in detail in previ- 
ous papers 0, 0] and will be discussed here only briefly 
to point out recent improvements. Our code evolves the 
full Einstein field equations coupled to relativistic hydro- 
dynamics in 3+1 dimensions. We have recently gener- 
alized this code using the Cartoon methods of [23, 
so that it can perform 2+1 simulations in axisymmetry 
in the same coordinate system. In order to improve its 
behavior near the intersection of the excision zone and 
the symmetry axis, we add a small amount of Kreiss- 
Oliger dissipation [53|, |5J, |55j to the evolution equation 
for the extrinsic curvature Aij. A further description of 
our axisymmetry algorithms, together with axisymmetry 
code tests, will be presented in a forthcoming paper |44j . 
in which the effects of viscosity on differentially rotating 
binary neutron star remnants are studied. 

We evolve the field evolution equations using the BSSN 
formulation H E3]. In the BSSN system, one decomposes 
the 3-metric as 7y = e 4 ^7y and the extrinsic curvature as 
Kij = e A ^(Aij + jijK/3), and one promotes the confor- 
mal connection coefficients f l = to independent 

variables. One then uses the ADM equations to write 
evolution equations for the new set of fundamental vari- 
ables: 7y, (f>, A^, K, and T\ On each time slice, these 



where AT is the timestep, chi — 0.1, ch2 = 0.5, and 
c H3 = 1 [For the complete right hand sides of Eqs. ©- 
©, see |23, equations (12), (11), and (14).] Modifica- 
tions similar to those in Eqs. 10 and Q were suggested 
in [60j, while a modification similar to Eq. JSJ has re- 
cently been used in for doing excision in the ADM 
formulation for pure vacuum spacetimes. Eq. J^J) intro- 
duces a diffusive term into the evolution of <j>. Eq. © 
introduces a nonlinear damping term into the evolution 
of A^. We find that modification (JSJ) has the largest 
impact on accuracy. 

Of crucial importance for the stability of our code are 
our constraint additions to the T 1 evolution equation. As 
shown in p4|. our equation for 9 t f l has the terms 



>),v ^r+\, F' r.j 



(9) 



Looking, for example, at the x-component of this equa- 
tion, 



(10) 



we see that if /3 J j > or f3 x ^ x < 0, then dtT x contains a 
term tending to produce exponential growth. We lessen 
the possibility of an instability caused by these terms by 
using J21 to replace (fTUfl with 



d t T x 



[Fj + \ A \Pj\] (-T^fcJ-^l^-lf* 



[d\ x + Ab |/3 x , 
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A B |/3%|r* 
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and similarly for T y and T z . Note that the "exponential" 
terms in the above equation (i.e. the terms proportional 
to r x ) are now guaranteed to be exponential decay terms. 
We find good results with A^ = 2/3 and As = 3/4. 

Alcubierre et al. [62] find improved behavior when they 
enforce the constraint T = 0. Yo et al. 0] found it useful 
to enforce T = T> = 0. We instead apply the reasoning 
above to modify the evolution equations for 7^ and Ay. 
Thus, in the equation for j xx , we find the terms 

A*.- (-!*, + (12) 

which we replace by 

2 1 ■ n 2 11 

dtTxa = g [-P 3 ,j + \P J G xx - -Ac ,j\ ixx 

+2 [/3% + A D \(3\ x \] G xx - 2X D \0 X , X \ % x 
+ ••• , (13) 

where G xx is the value of j xx as computed from the five 
other independent components of 7y, assuming I? = 0. 
We perform the same substitution for ^ yy and ^ zz . We 
use Ac = 2/3 and Ad = 1/10. In a similar fashion, we 
modify the evolution of A xx , A yy , and A zz from 

d t A xx = ■■■ + (-\p 3 ,j + W\x + aK^j A xx (14) 

to 

2 2 - 

d t A xa; = - [-/3 J j + Ac \f3 3 ,i|] ff xa; - -Ac |/? J j I A xa; 

+2 [/3% + A D |/3%|] H ax - 2A D |r ,x| A ra 

+ [aK + \ E \aK\] H xx - 2\ E \aK\ A xx 

+ ••• , (15) 

and similarly for the other two components. Here Xe = 
0.1, Ac and Ad are the same as above, and H xx is the 
value of A xx computed from the five other independent 
components of Aij assuming T = 0. 

We take spatial derivatives in a centered way — we do 
not use causal differencing. The only exception, as sug- 
gested by |ri|. is in the advection terms along the shift 
f3 1 di 1 for which we use the second-order upwind differenc- 
ing described in |6^| . 

Our hydrodynamics scheme uses van-Leer type advec- 
tion and artificial viscosity shock handling [24[. It is 
known that such schemes can be inaccurate for ultra- 
relativistic flows [64|]. We monitor the Lorentz factors of 
our fluids, and find that they never exceed «2, which is 
around the upper limit for accurate evolutions with a van 
Leer code. In addition, most of our runs do not involve 
strong shocks. We thus believe that our hydrodynamics 
scheme is adequate for the present purposes, although wc 
eventually may have to improve it. Our hydrodynamics 
scheme employs the "no atmosphere" approach |24j . so 
that the density at any point on our grid is allowed to 
fall to zero. It is important that we are able to dispense 



with an artificial atmosphere. If we could not, then in 
situations where all the matter in the problem falls into 
the black hole, the hole would continue to accrete atmo- 
sphere indefinitely, and its mass would continue to grow 
unphysically. 

The boundary conditions we apply at the edge of the 
excision zone are described in detail in |l2^ . They consist 
of taking the time derivatives of quantities at the excision 
boundary from the time derivatives of these quantities at 
adjacent points. We use spherical excision regions inside 
the apparent horizon throughout (see |65j regarding the 
superiority of spherical to cubic excision regions). We 
have tried several boundary conditions for the matter 
variables, and have found that our results arc insensitive 
to the choice, as they should be. In the runs described 
below, we simply set the matter variables equal to zero 
when they hit the excision zone, thus making the excision 
boundary a perfect one-way membrane. 

The lapse and shift must be chosen in such a way that 
the total system of evolution equations is stable. It is 
also desirable that the gauge conditions are chosen so 
that, as the system settles into equilibrium, it appears 
stationary in the adopted coordinates. We have experi- 
mented with several choices for the lapse a and shift f3 z , 
and we have found that driver conditions using the sec- 
ond time derivatives of a and f3 l provide the most stable 
evolutions. Following the suggestion of Alcubbierre et 
al [13, we have had great success with the hyperbolic 
shift driver condition: 

did 1 = b^adtT 1 - b 2 d t f3 l ) , (16) 

with 61 = 0.75 and b 2 = 0.27A/- 1 (c.f. 11]). One can 
create a hyperbolic lapse condition by introducing two 
coupled first-order equations and a new function A 

dta = aA 

d t A = -ai(ad t K + a 2 d t a) , (17) 

with ai = 0.75 and a 2 = 0.27M -1 . The a in front of 
A in the first equation is a "safety" feature, to prevent 
the lapse from dropping to zero. With this safety fea- 
ture, we find that the lapse levels off at finite positive 
values everywhere on and outside the excision zone for 
all our runs, thereby maintaining a "horizon penetrat- 
ing" (a > 0) time coordinate. However, at late times 
(i ~ 200M), we find that the asymptotic values of some 
of our variables (e.g. ^ xx ) begin to drift, increasing lin- 
early with time. This drift is also present when harmonic 
slicing, another slicing with a hyperbolic character |66| . 
is adopted. Apparently, Eci. (|17|) does not sufficiently re- 
strict the coordinate system's evolution. We remove the 
drift by adding a third term to Ea. l|17|) proportional to 
K— -Kdrivo where -ftTdrive is some reasonable positive func- 
tion. In this way, the value of K itself, and not just its 
time derivative, is "driven" . We shall refer to this slicing 
as our "hyperbolic lapse" . The complete slicing condition 
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presented in |l2j. The final integrals are 



dta = aA 

d t A = -a x {ad t K (18) 
+a 2 {dta + e- 4 *a(K~K drivc )}) . 

Here the e _4 ^a factor is chosen so that the new term is 
small in the strong-field region, where (|17|l works well, 
but becomes comparable to the other terms in the outer 
portions of the grid, where it successfully removes the 
drift. 

We have tried several forms for i^drivc- The simplest, 
and usually adequate choice, is zero. This drives K to 
zero (maximal slicing) and usually causes a very slow 
downward drift in the lapse near the horizon. For many 
astrophysical applications, where we only need to evolve 
for several hundred M , this is usually unimportant. How- 
ever, the effect can be removed by a better choice of 
-f^drivc- One possibility is -Kmit, the value of K at the 
time excision is introduced. Another choice is Kks, a 
function whose form is inspired by the Kerr-Schild rep- 
resentation of a Kerr black hole [c.f. Eq. (36) of [T^]. 



H 

V 



2a 3 (l + H)l i H <i + 2aHl\ i (19) 



1 



1) 



2 

/37(2a 2 J?) 



Note that when we choose this functional form for Kks, 
the lapse and shift typically are not the same as the Kerr- 
Schild a and [3 l . 

For K = -Kdrivc, we apply our usual excision boundary 
conditions on a. Otherwise, there are no spatial deriva- 
tives in Eq. (|18fl . and no explicit inner boundary condi- 
tion is needed. In some cases, however, we have found 
more accurate results when we hold the values of the 
lapse on the excision zone fixed in time (the "frozen" 
inner boundary condition). 



M = 



16tt 



d 3 3 
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16ttp + AijA* - -K 2 
-t^ k f jik + (1 - e+)R 



-—<b {r-8D^)dSi 



(22) 



8tt 



e b0 ( Ah + ^xW k K 



--x'A in d k f n + 8iicx>ak 
e^x^A'kdSi. . 



d 6 x 



(23) 



on 



We choose the inner surface 9f2 to be a sphere with a co- 
ordinate radius about twice that of the excision bound- 
ary. This puts dfl slightly outside the apparent horizon 
in the simulations reported below. 

In our pre-excision code, f2 is chosen to cover the entire 
numerical grid, and there is no surface integral contribu- 
tion. The rest mass Mq cannot be used as a diagnos- 
tic because it is conserved identically in our pre-excision 
code. Our pre-excision code also conserves J identically 
in axisymmetry |44j . With excision, Mo is not expected 
to be conserved in fl, since matter falls into the excision 
region. When evolving with excision, J is not identi- 
cally conserved, even in axisymmetry, and thus serves as 
a code check together with M. 

Once a black hole is present, we detect it by using an 
apparent horizon finder (see j6g for details). As the sys- 
tem approaches stationarity, the apparent horizon will 
approach the event horizon. We estimate the size the 
horizon in our coordinate system by the radius tah con- 
structed from the I = 0, m = moment of the horizon 
surface. From the surface area of the apparent horizon, 
we compute the irreducible mass Mj rr defined by 



(24) 



III. DIAGNOSTICS 

Our most important diagnostics are the conserved 
mass M and angular momentum J. These are both de- 
fined by surface integrals at infinity [67| |: 



M = -1- / ^ll tm i Jn {lmn. ] -l ] n. m )d 2 S l (m 

Ibir J r=00 



J - — £ ' 



x 3 K™d ,L 



(21) 



We also compute the proper circumference of the hori- 
zon in the equatorial (xy) plane, which we call C eq , and 
we compute the proper circumference in the meridional 
(xz) plane, which we call C po \. For static nonrotating 
black holes, C cq = C po \ — AwM. For stationary rotating 
black holes, one can compute C eq and C po i from the Kerr 
metric in Boyer-Lindquist coordinates to be 



C cq = AttM 



(25) 



,tt/2 i 

Cpoi = 4M / d0d2 + 2y/l -q 2 + q 2 sin 2 6 ,(26) 
Jo 



We measure M and J by applying Gauss' Law to obtain a 
surface integral over an inner surface, <9f2, (which encloses 
the singularity) plus a volume integral over the space 
outside this surface, f2. Details of this calculation are 



where q = J/M 2 is the spin parameter of the black hole. 
The ratio C po i/C eq varies from 1 for q = to 0.6 for 
q = 1. For the black holes in our simulations, we infer 
the horizon mass Mah from C eq and Eq. 1)25(1. We infer 
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FIG. 1: The evolution of the mass M, angular momentum 
J, and lapse variation Aa for the evolution of an a/M = 0.4 
black hole in Kerr-Schild coordinates. We use a 30 2 grid to 
cover the meridional plane.. 

the horizon angular momentum Jah from C po i/C eq and 
Eq. {2H>, together with Mar- 

Finally, we find the ergosurface of the black hole. The 
ergosphere is defined in the stationary limit, in which 
case -Q^ is a Killing vector, and the ergosurface is defined 
as the surface where 500 = §1 • = 0, with g 00 > inside 
and 500 < outside. 

As in [H 0, we gauge the degree to which a field 
/ reaches stationarity by monitoring Af(t), defined to 
be the L2 norm of f(t) - f(t- AT), where AT is the 
timestep. We compute the L2 norm of a gridfuntion g by 
summing over every gridpoint i: 



^%) = JE^ (27) 



IV. TESTS 
A. Field Code Test: Vacuum Black Holes 

In a previous paper [l^, we used our code to evolve 
isolated, stationary black hole spacetimes in Kerr-Schild 
coordinates. These coordinates have the advantages of 
being horizon-penetrating (a ^ at the horizon) and 
providing a manifestly stationary metric. We were able 
to evolve both stationary and rotating black holes for ar- 
bitrarily long times. We succeeded in doing this both 
when evolving only one octant of the space and when 
evolving the full space without any symmetry assump- 
tions. These evolutions were done in three dimensions 
using a different set of gauge conditions from those uti- 
lized in this paper. In Figure ^ we show the evolution 



of a a/M = J/M 2 = 0.4 Kerr black hole in Kerr-Schild 
coordinates using our 2D axisymmetry code and our hy- 
perbolic gauge conditions. For this case, we use a frozen 
inner boundary condition on a, and turn off the third 
term in i|18|) . (Using i^drivc = ^init gives similar re- 
sults.) We use a grid spacing of AX = AM, with outer 
boundaries at 12M and an excision zone at a coordinate 
radius of 1.5M, as was used by [l2(. The event hori- 
zon for a/M — 0.4 is located at r oq = 1.917M in these 
coordinates. 

As a second test, adapted from we evolve a 

Schwarzschild black hole in initially isotropic coordinates. 
Choosing a = 1 and (3 % = at t — 0, the initial metric is 

ds 2 = -dt 2 + (l + ^ (dx 2 + dy 2 + dz 2 ) , (28) 

where r = \/ x 2 + y 2 + z 2 . The event horizon is located 
at r = 0.5M in these coordinates. Physically, this black 
hole is stationary, but it does not appear stationary in 
the coordinates generated by Eqns. (|16|) and (|18l) starting 
with the initial lapse and shift cited above. By evolving 
this spacetime, we check that our excision code can work 
with coordinates other than stationary Kerr-Schild. We 
also check the ability of our gauge conditions to "find" co- 
ordinate systems which make the metric manifestly sta- 
tionary. We allow the lapse to drop, so we do not freeze 
the lapse at the excision zone, but employ equation 118|) 
everywhere. We use i^drivc = K mi t = 0, since Kks is 
singular for our value of a at t = [see equation i|19|) ]. 

In Figure |3 we plot the results for a run in axisymme- 
try with outer boundaries at 12M, an excision radius of 
0.36Af, and a grid of 128 2 to cover the meridional (xz) 
plane. Also shown are scaled results for a 64 2 run to 
demonstrate convergence. We also performed a run on a 
256 2 grid with the same resolution as the 128 2 run but 
with the outer boundaries at 24M. From the figure, we 
see that the error can be controlled by the grid resolution 
and the location of the outer boundaries. We see that the 
surface area of the apparent horizon (i.e. M- m ) remains 
nearly constant while the coordinates adjust to create a 
stationary system. This indicates that the apparent hori- 
zon is following the event horizon well. The coordinate 
adjustment is reflected in the initial increase in the coor- 
dinate radius of the horizon and in the drop of the lapse. 
Note that the lapse settles quickly, and that it remains 
positive everywhere outside and at the excision zone. To 
check that the black hole remains a Schwarzschild black 
hole, we monitor C oq and C po \ and find that they both 
remain equal to 47rM to within one percent. 

B. Hydro Code Test: Relativistic Bondi Flow 

Next, we test our hydrodynamics code by solving an 
accretion problem that has an exact solution. In a pre- 
vious paper |24| . we confirmed our code's ability to ac- 
curately simulate shocks, spherical dust collapse, nonro- 
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FIG. 2: The evolution of a nonrotating black hole in our 
hyperbolic gauges, starting in isotropic coordinates with a = 
1, /3 Z — 0. On top, we show the deviations of the ADM 
mass Madm and the irreducible mass Mi rr from their initial 
value: SM = (M — Mi) /Mi. SM is shown for runs with outer 
boundaries at 12Mi using a 128 2 grid and using a 64 2 grid, 
to demonstrate convergence. We also show a run with outer 
boundaries at 24M; using a 256 2 grid to determine the effect 
of the outer boundary. Below, we show the time evolution on 
the 128 2 grid of the apparent horizon coordinate radius tah 
and the maximum values of a and j X x on the grid. 

tating and rotating polytropes, and binary polytropes. 
Now we test its ability to maintain stationary, adiabatic, 
spherically symmetric accretion onto a Schwarzschild 
black hole, in accord with the relativistic Bondi accretion 
solution for T = 1.5 69] . Following the suggestion of [7pj . 
we write the metric in Kerr-Schild (ingoing Eddington- 
Finklestein) coordinates; in this way, all the variables 
are well behaved at the horizon. We begin by holding 
the field variables fixed in order to prevent the black hole 
from growing due to accretion. 

We evolve this system twice, once using a 64 2 grid in 
2+1 and once using a 64 3 grid in 3+1. We place outer 
boundaries at Y1M and an excision zone at a coordinate 
(areal) radius of 1.5M. At t = 0, we set the density 
and velocity profiles according to the exact solution for 
r = 1.5 and an accretion rate dM /dt\ &cc — 0.0031, with a 
sonic radius at 10 5 M. This accretion rate is maintained 
through the evolution by fixing the hydrodynamic vari- 
ables on the outer boundaries at their exact steady-state 
values. In Figure El we plot A(po), as defined in Sec- 
tion IIIII and also the values of po at selected points in 
the accretion flow. For A(po), we reach machine pre- 
cision after less than lOOAf , making further integration 
unnecessary (The velocity fields have also frozen 

near their initial values by this time.) 

When we allow the fields to evolve, we see the ir- 
reducible mass of the hole grow at a rate dM- m /dt w 



FIG. 3: The settling of the rest-mass density to steady-state, 
starting from the analytic value. The change per timestep 
quickly drops to the machine level. On top, we plot Apo for 
both the 64 2 2D run and the 64 3 3D run. Below, we show 
the time evolution of po at three points on the diagonal line 
x — y — z in the 3D run, each normalized to its initial value. 
pi corresponds to po measured at r = V3x = 2M, p2 to po at 
r = 6M, and p3 to po at r — 10M. 



0.9dM/di|acc- This error is consistent with the errors in 
our irreducible mass found at this numerical resolution, 
even in the absence of accreting matter. 



C. Oppenheimer-Snyder Collapse 

Next, we simulate the Oppenheimer-Snyder collapse 
of a homogeneous spherical ball of dust to a black hole. 
The behavior of this system is known in several coordi- 
nate systems [Z^, El ■ We use a 160 2 grid with outer 
boundaries at 1AM. At t = 0, the dust is at rest and 
has an areal radius of 3M. We start in an isotropic co- 
ordinate system, in which 7^ = 8ij. Our initial a and 
(3 l are set by enforcing maximal slicing and the minimal 
distortion gauge condition, respectively (see [z3)- Since 
the ball has no pressure support, it immediately begins 
to collapse. During the first phase of this collapse, there 
are no trapped regions and no singularities, so we evolve 
the entire grid without excision. Our code checks du ring 
this part of the evolution are well satisfied; see [2J, [44| • 
For gauge conditions during this no-excision phase, we 
use our hyperbolic lapse and shift drivers. We evolve 
in this way from t = to t = HAf, at which point 
our no-excision code crashes due to its inability to re- 
solve the central region ("grid stretching"). An apparent 
horizon appears at t = 9M at a coordinate radius of 
^ah = 0.96M with an irreducible mass of M- m = 1.02M. 
We next repeat the evolution from t = 10M with our 
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FIG. 4: ADM mass, horizon diagnostics, and Aa for the 
collapse of a homogeneous sphere of dust to a Schwarzschild 
black hole. Collapse begins at t = and black hole excision 
occurs at t = 10M. 



FIG. 6: ADM mass, horizon diagnostics, and Aa for the col- 
lapse of a nonrotating, unstable n = 1 polytrope from appar- 
ent horizon formation at t/M = 27 through final stationarity. 
The code is axisymmetric and uses a 128 2 grid. 
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settle quickly. As we continue the evolution, the remain- 
ing exterior rest mass falls into the excision zone over the 
course of the next 100M, and we are left with a vacuum 
spacetime. We evolve for 400M, by which time the sys- 
tem has long since settled to a Schwarzschild black hole. 
Oppenhcimer-Snyder collapse does have an analytic solu- 
tion in Friedmann coordinates, but not in the coordinates 
we are using, which are defined by our gauge conditions 
(|16f) and i|18|l together with the boundary conditions on 
a and [3 l at r ex . Therefore we check the accuracy of our 
evolution using global invariants. In Figure 01 we show 
our mass diagnostics for the post-excision run, which con- 
firm that the end product is a Schwarzschild black hole, 
and we plot Aa as proof of stationarity. In Figure [5] we 
plot the magnitude of the constraint violations as func- 
tions of time. These show that the error is not growing 
during the long stationary evolution. 



FIG. 5: Violation of the Hamiltonian Tt, momentum M % , 
and Gamma Q l constraints as a function of time for the 
collapse depicted in Fig. 21 We plot the un-normalized L2 
norms, where we use the shorthand L2(M i ) 2 = L2{M X ) 2 + 
L2{M y ) 2 + L2{M Z ) 2 and 

L2{g i f = L2{g x f + L2{g v ) 2 + L2(g z f. 



excision algorithm and an excision boundary at radius 
Tex = 0-7 M. At this point, only 1.2% of the rest mass 
is outside the horizon, but the spacetime in our coordi- 
nates is still changing. We continue to evolve with our 
hyperbolic gauges, and we allow a to drop at the excision 
boundary. In this example, using -Kdrivc = Kks is far su- 
perior to any other choice, since only then does the lapse 



D. Collapse of a TOV Star 

The previous example possessed spherical symmetry 
and no pressure. In our next test, we study the collapse 
of an unstable nonrotating, spherical polytrope, whose 
initial state is given by the solution to the TOV equa- 
tions |7^ . 

For initial data, we take a perfect fluid with equation 
of state P = npo 1+1 / n , with n = 1, and we choose our 
units such that K = 1 75] . In these units, the n = 1 
TOV sequence has a turning point at the critical central 
rest density p" lt = 0.32 where the ADM mass of the 
star is M max = 0.164. We choose to evolve a star with 
initial central rest density p c = 0.5 and ADM mass M = 
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FIG. 7: Same as for Figure |SJ but now the collapse is simu- 
lated on a 3D 64 3 grid. 

0.158. As this star is on the unstable branch of the n — 
1 sequence, it is unstable to radial oscillations and will 
collapse to a black hole. We evolve the first part of the 
collapse without excision using a 128 2 grid, with outer 
boundaries at 12. 7M and with our hyperbolic drivers. 
We evolve from t = to t = 28. 5M, locating an apparent 
horizon at t = 27M with radius r = 0.6M and irreducible 
mass Mj rr = 0.95A/. We begin an excision run from 
t = 27. 8M, at which point 4% of the rest mass is still 
outside the apparent horizon and 8% is outside of the 
excision zone. All of this matter falls into the excision 
zone by t — 31.6M. It should be emphasized that the 
spacetime in these coordinates is more dynamical than 
the above numbers might suggest: e.g. during the first 
10M of post-excision evolution, the maximum value of 
A^AijM 2 increases from 0.25 to 0.44. The system settles 
quickly thereafter, as we see by evolving an additional 
350M to 390Af . In FigureEl we show our diagnostics for 
this run. 

All the runs described above were carried out on two- 
dimensional axisymmetric grids. In Figure we show 
diagnostics for the same collapse in a three dimensional 
simulation, with a 64 3 grid and boundaries at [0, 12. 7M] 3 
(employing octant symmetry to evolve only the upper 
octant). The behavior of each quantity is similar to that 
in the 2D run. 



E. Collapse of a Rotating Star 

Gravitational collapse of astrophysically realistic stars 
will involve rotation. Even if the progenitor star rotates 
slowly, it will spin up as it collapses if it conserves angular 
momentum. It is therefore important to test our code by 
simulating the collapse of a rapidly rotating star. 



FIG. 8: Mass M and angular momentum J during the post- 
excision phase of the collapse of star A. We show results for 
axisymmetric runs carried out with a 80 2 , a 160 2 , and a 320 2 
grid. Both M and J are measured in two ways. The solid lines 
are quantities as measured by the integrals (I22li and (IL'MI . The 
dashed lines are obtained by measuring the geometry of the 
apparent horizon and comparing with the Kerr metric (Mah 
as inferred by C eq and Jah by C vo \/C eq ). For the 160 2 and 
320 2 runs, the two J measurements lie on top of one another. 



The star we adopt as initial data, labeled A, is de- 
scribed in Table [I] The initial data was obtained using 
the relativistic equilibrium code of [j^. Star A is a "hy- 
permassive star" with a mass M — 0.19, which is 20% 
higher than AT max , the maximum allowed mass of a non- 
rotating TOV star. Star A is able to maintain this mass 
because of the added support against gravity provided 
by (differential) rotation. The star has J/M 2 = 0.57, so 
that the eventual Kerr hole will have appreciable spin, 
assuming all of the mass and angular momentum is cap- 
tured by the hole. Even prior to collapse, the effects of 
angular momentum on the star are significant, as we can 
see by noting that the radius of the star on the rota- 
tion z-axis (the polar radius) is only 70% of the radius of 
the star in the equatorial plane (the equatorial radius). 
Star A has a differential rotation profile (see next sec- 
tion), so there are no turning-point theorems which can 
be applied to determine the stability of this star, but we 
find numerically that it is unstable to collapse. Pertur- 
bations due to numerical (roundoff) error are sufficient 
to trigger the collapse, but the onset timescale for col- 
lapse is not independent of resolution. In order to do 
convergence studies, we deplete a small percentage (4%) 
of the initial pressure, so that the initial perturbation 
is resolution- independent. This perturbation is so small 
that re-solving the constraint equations at t = makes 
little difference. 

We carry out the entire evolution, before and after exci- 
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sion, in the hyperbolic gauges. (The choice of i^drive has a 
negligible effect on the evolution in this application.) We 
perform the same evolution on a 80 2 grid, a 160 2 grid, 
and a 320 2 grid. On the 320 2 grid, a horizon appears 
in the pre-excision run at t — 44M, with instantaneous 
radius tah = 0.5M and mass M- lrr — 0.77M, which are 
growing rapidly. We excise at time t = 45. 5M and radius 
r ox = 0.43M, so that 22% of the total rest mass is still 
outside the excision region, and 12% is still outside the 
apparent horizon (which now has radius tah = 0.73A/). 
This matter quickly falls into the hole, and, after evolving 
for 6M with excision, the exterior spacetime becomes a 
vacuum. In Figure [3J we check the ability of our code to 
conserve mass and angular momentum during this phase 
of the evolution. The mass is well conserved on all three 
grids, but the angular momentum slowly decreases with 
time. Increasing resolution reduces this loss of J. The 
violations of the constraint equations also converge to 
zero as resolution is increased. We can evolve stably for 
t ^S> 100M, but the loss of angular momentum is too great 
past this point for the evolution to be reliable unless the 
grid exceeds 320 2 . 

Figure [S] suggests that the angular momentum loss can 
be controlled by increasing resolution. Moreover, we have 
already shown that our code can conserve J for an ar- 
bitrarily long time while evolving a Kerr black hole in 
Kerr-Schild coordinates (see Fig. Q] and 01 )• Given this 
fact, we could eliminate the J loss by transforming to 
Kerr-Schild coordinates when we introduce excision. Al- 
ternatively, we might carry out the entire evolution in 
Kerr-Schild-like coordinates. (This would require devel- 
oping gauge conditions which would force the coordinate- 
system to maintain its Kerr-Schild-like character as the 
system evolves.) We are currently investigating these 
possibilities. In the meantime, we can already evolve such 
matter-black hole systems long enough to tackle several 
interesting problems. 



V. APPLICATION: THE COLLAPSE OF 
RAPIDLY ROTATING STARS 

Tracking the collapse of rapidly rotating stars is one 
of the most important applications of numerical general 
relativity. Such simulations determine the fate of col- 
lapse and provide a test of the cosmic censorship conjec- 
ture [73 • If the star collapses to a stationary black hole, 
the "no-hair" theorems require that it settle down to a 
Kerr black hole. In the Kerr spacetime, the singularity is 
covered by an event horizon only if q = J/M 2 < 1; other- 
wise the singularity is naked. Rotating stars, on the other 
hand, are not so restricted, and sufficiently rapidly rotat- 
ing stars will have q > 1. When these stars collapse, it 
thus seems conceivable that they could form naked singu- 
larities. Alternatively, if the cosmic censorship hypoth- 
esis [?3 is true, then the collapse of the whole system 
must somehow be averted. This can happen if the star 
loses angular momentum as it collapses, either by gravi- 



TABLE I: Equilibrium Star Configurations 
(n = 1, Mo = 0.2). 



Star 


M a 


R b 

TLcq 


Rc c 


1 d 


T/\W\ e 




Fate h 


A 


0.19 


0.6 


0.8 


0.57 


0.10 


0.29 0.70 


BHND 


B 


0.19 


1.2 


1.4 


0.91 


0.18 


0.38 0.50 


BHND 


C 


0.19 


1.6 


1.8 


1.18 


0.23 


0.40 0.39 


NBH 



a ADM mass 

b coordinate equatorial radius 
c areal radius at the equator 
d q = J/M 2 

e ratio of kinetic to gravitational potential energy 

* ratio of central to equatorial angular velocity 

9 ratio of polar to equatorial coordinate radius 

h BHND = black hole, no disk; NBH = no black hole 

tational wave emission or by shedding matter with high 
specific angular momentum, so that the final black hole 
has q < 1. A naked singularity can also be averted if the 
collapse of a q > 1 star is always halted by centrifugal 
forces, so there will be no black hole and no singularity 
at all. Nakamura [j^ has pointed out that a centrifugal 
barrier could protect cosmic censorship in this way. As- 
suming no mass or angular momentum are shed during 
the collapse, the radius Rb at which the centrifugal force 
balances the gravitational force will be 



R 2 M 2 R? b ' 

so that 

R b ~ Mq 2 . (30) 

Nakamura argues that, if q < 1 (i.e., the star is subKerr), 
the star will already be inside a black hole before rota- 
tion can halt the collapse. For q > 1 (i.e., the star is 
supraKerr) , the collapse will be halted at a radius larger 
than M, and no black hole forms. 

Shapiro and Teukolsky 54] have studied the collapse 
in full general relativity of axisymmetric tori consisting 
of collisionless matter, and have found that black holes 
form only from subKerr initial configurations. The first 
numerical simulations of the collapse of rotating relativis- 
tic fluid stars were carried out in axisymmetry by Naka- 
mura |73 and Nakamura and Sato [8(J • They found that 
a black hole forms only when a subKerr star collapses. 
(For stars with q within 5% of the critical value, Naka- 
mura |79j could not determine the final fate and could not 
exclude the possibility of a naked singularity.) Stark and 
Piran |8lj also performed simulations which showed q~l 
to be the critical point of demarcation between collapse 
and bounce. Shibata |53 performed a detailed study of 
the collapse and bounce of subKerr stars in axisymme- 
try. These hydrodynamic studies did not (and sometimes 
could not) study in detail the fate of the matter in the 
outer layers of the star when a black hole forms. More 
recently, Shibata [2^ has studied the collapse to black 
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FIG. 9: Snapshots of the rest-density contours and the velocity field (v x , 

star B to a black hole. The contour lines are drawn for p = 10 - ' ' 2 p M ax f or j = o, 1, .., 12. Prior to excision, po' ax is set 

equal to the instantaneous maximum value of po. Afterwards, it is held at the maximum of po at the time of excision. Vectors 
indicate the local velocity field, v*. The thick curve in the last three frames marks the apparent horizon. On the last frame, 
the exterior spacetime is nearly a vacuum. 



holes of uniformly rotating polytropes spinning at the 
mass-shedding limit. He finds that, for polytropic indi- 
cies 2/3 < n < 2, the star collapses to a Kerr black hole 
with no appreciable disk. By using high resolution, he 
is able to follow the system for At ~ 20M after an ap- 
parent horizon is first located. This time approaches the 
limit of reliable evolution without excision, but in this 
case it is long enough to see all the matter fall into the 
hole. By contrast, Shibata and Shapiro 2JJ considered 
the collapse of an n — 3 polytrope spinning uniformly at 
the mass-shedding limit. Such a configuration is nearly 
Newtonian (R cq — 620M) at the onset of collapse, and 
it forms an appreciable disk (Md/M w 0.1) around the 
final black hole. While the final disk mass can be esti- 
mated from the angular momentum distribution of the 
outermost regions (see also H2)> an d a l so by extrapolat- 
ing the growth of the black hole horizon to late times, it is 
not possible to follow the final relaxation to a stationary 
state without excision or to probe for nonaxisymmetric 
instabilities that may arise in the ambient disk |47j 

Our excision code should be well suited to finding the 
final state of any rapidly rotating stellar collapse — not 
only for determining whether or not a black hole forms, 
but also for determining how much rest mass escapes 
collapse if one does form. To explore this capability, 
we take differentially rotating polytropes as our initial 



data, so that we can study both subKerr and supraKerr 
cases. Differential rotation is naturally produced in su- 
pernova core collapse |8^|, accretion induced collapse of 
white dwarfs to neutron stars |8^ . and binary neutron 
star coalescence 0, E2> ■ Our adopted rotation law 
is 

u'u^^ R 2 cq A 2 {n c ^n) , (31) 

where f2 is the angular velocity of the fluid, f2 c is the 
value of f2 on the rotation axis, R cq is the equatorial co- 
ordinate radius. The parameter A measures the degree of 
differential rotation and is chosen to be unity for all cases 
below, so that the centers of our stars rotate about three 
times faster than their equators. We take the z-axis to 
be the rotation axis, and define the cylindrical coordinate 
radius w = y/ ' x 2 + y 2 . In the Newtonian limit, Eq. 131(1 
reduces to the so-called "j-constant" law [86j 



We choose a polytropic index n = 1, and take our ini- 
tial stars to be sufficiently compact so that the collapse 
does not span a large dynamic range. Accordingly, we 
are able to use a single, modest grid for each run. As 
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in [23 . l8l| , we induce collapse by depleting the initial 
pressure by a factor: P — > /pP. Below, we show results 
for /p = 0.01. While this form of artificially-induced col- 
lapse does not correspond to any realistic astrophysical 
scenario, there are several situations in which an "effec- 
tive" pressure depletion does occur. For example, the 
collapse of the core of a massive star which produces a 
supernova is brought about by the removal of pressure 
support both from photo-dissociation of iron-nickel nu- 
clei and the neutronization of the core (de-leptonization) . 
Phase transitions in neutron stars, such as a transition 
to quark matter, or rapid de-leptonization via neutrino 
cooling, could also have the effect of inducing pressure 
depletion. We choose fp = 0.01 to make pressure forces 
unimportant in comparison with centrifugal forces and 
gravity. After depleting pressure from the star, we re- 
solve the constraint equations to produce valid initial 
data. This process of depleting pressure and re-solving 
the constraints causes M and J to drop by a few percent, 
while J/M 2 changes by one percent or less. 

Table I lists the equilibrium stars used to construct our 
initial data. These initial data were generated using the 
code of . Each star has the same rest mass M — 0.2, 
so our stars are members of a sequence uniquely defined 
by n = 1, A = 1, Mo = 0.2. This sequence crosses q = 1 
at one point, between our second and third stars, stars B 
and C. We expect to find a qualitative difference in the 
behavior of stars B and C. 

Star A is exactly the star studied in the previous sec- 
tion. It is dynamically unstable and collapses without 
pressure depletion to a Kerr black hole with no disk. Not 
surprisingly, this is also found to be the behavior when 
pressure is depleted. We will concentrate below on stars 
B and C. We begin with simulations in axisymmetry and 
then discuss simulations in full three dimensions. 
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FIG. 10: Diagnostics for the collapse of star B. Above, we 
show the evolution of M and J calculated from integrations 
of the exterior spacetime and from measurements of the ge- 
ometry of the horizon. Below, we plot the total rest mass on 
the grid, normalized to its initial value. Rest mass is con- 
served prior to excision. At t — 30M, we excise a region from 
the middle of the grid. This cuts out the matter inside this 
region, which accounts for 80% of the total rest mass. Over 
the next 20M, the remaining rest mass falls into the excision 
zone, leaving a vacuum being evolved in the outside region. 



tern has, however, been entirely determined well before 
this time. 



B. SupraKerr Collapse 



A. SubKerr Collapse 

Star B has J/M 2 = 0.9, so it is subKerr. Its collapse in 
axisymmetry is shown in Fig. EI We evolve on a 300 2 grid 
with outer boundaries at 14M. At t = 28AM, we locate 
an apparent horizon with tah = 0.62M, M- ltI = 0.72M. 
We excise at t = 29M, at which time M irr = 0.74M, 
22% of the rest mass is outside our excision zone, and 
15% is outside the apparent horizon. The horizon cir- 
cumferences at this time are in the ratio C po i/C cq = 0.76, 
which, if this were a stationary Kerr horizon, would cor- 
respond to q = 0.92. We continue evolving with an exci- 
sion boundary at radius r cx = 0.08. All of the matter falls 
into the hole within 20Af after excision is introduced. We 
evolve for an additional 20M after this. We find no signs 
of numerical instability. Mass conservation is excellent 
(the amount lost due to gravitational radiation is below 
0.1%), but the gradual loss of angular momentum noted 
in Section llV El is present, as can be seen in FigurelTUl We 
stop evolving when the total angular momentum drops 
below 80% of its initial value. The final state of the sys- 



Star C has J/M 2 = 1.2. We remove the star's pressure 
support and evolve. In Figure lTTl we show the results of a 
400 2 axisymmetric run with boundaries at 13 M. With its 
pressure support removed, the star immediately flattens 
along the z-axis and moves inward in w. This inward 
motion toward the axis is halted by centrifugal forces. 
As seen in the upper right panel of Fig. 1111 the inner 
region of the star stops collapsing before the outer re- 
gion, so a strong shock is formed. The star then expands 
into a torus whose radius oscillates with a period close to 
the initial central rotation period. We show the effects 
of this oscillation on the maximum rest density and the 
minimum lapse in Fig. 1131 We follow the torus for three 
oscillations during which time all our constraints are sat- 
isfied to better than 10%. The angular momentum J 
is conserved identically by our no-excision axisymmetric 
code, but we do find that the ADM mass M decreases 
gradually with time. This decrease cannot be accounted 
for by the small flux of rest mass and gravity waves out 
of the computational domain; the loss therefore repre- 
sents numerical error. We stop our evolution after three 
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FIG. 11: Snapshots of the rest density contours and the velocity field (v x , v z ) in the meridional plane during the axisymmetric 
collapse of star C to a torus. The contours are set as in Figure [5] Some velocity arrows appear outside the contours because 
the density there is very small but nonzero. Time is normalized to the initial central rotation period of the star, P ro t.c = 98M. 



oscillation periods because M has decreased by ~15%. 
To check that the evolution is qualitatively correct, we 
performed the same run on a 200 2 grid and found that 
the collapse, torus-formation, and oscillation of the star 
are very similar at this resolution. 

The torus formed in the above simulation could be 
subject to various non-axisymmetric instabilities. If the 
rotating torus fragments, the system may produce a large 
gravitational wave signal ("splash radiation" |E3)- ^ 
is therefore necessary to perform the above simulation 
in 3+1 dimensions. We perform this simulation using a 



280 x 140 x 200 grid, with boundaries at [-13M, 13M] x 
[0, 13M] 2 , where we use equatorial and 7r-symmetry. The 
results are shown in Fig. [21 The collapse, flattening, and 
formation of the torus occur as in the 2D runs. Then the 
torus quickly fragments into four clumps symmetrically 
located about the origin, roughly 90° apart. As these 
clumps collapse, they ultimately become too small to be 
evolved accurately on our grid. We conserve M and J 
to better than 10% throughout the integration shown, 
and we terminate the calculation when our errors exceed 
these bounds. To check this result, we have performed 
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FIG. 12: Snapshots of the rest-density contour lines for po and the velocity field (v x ,v v ) in the equatorial plane for the 3D 
collapse, bounce, and fragmentation of star C. The contours and time normalization are set using the same rule as in Figures |5] 
and llll Note that the origin of the system is now shifted to the middle of the rr-axis in this plot. 




FIG. 13: The maximum value of po and the minimum value 
of a during the evolutions of star C on different grids, plot- 
ted as a function of the initial central rotation period P ro t,c- 
The two 2D runs are qualitatively similar. The 3D run be- 
haves similarly to the 2D runs for about the first 0.5P r ot,c- 
Thereafter a nonaxisymmetric instability develops, and the 
collapsed star fragments. 



the same run on 140 x 70 x 180 and 100 x 50 x 100 grids. In 
each case, the torus fragments into four pieces 90° apart. 
In Fig. 1131 we compare the behavior of the maximum of 



Po and the minimum of a for the evolution of star C in 
3D to their behavior in 2D on 200 2 and 400 2 grids. 

It has been pointed out by Truelove et al. [§3 that 
spurious fragmentation may occur in a numerical simu- 
lation if the Jeans length is not well resolved. The Jeans 
length is given by 



(33) 




where p is the density (Mass-energy density and rest- 
mass density are nearly equal.) and c s = is the 

sound speed. We can get a lower bound on Aj by ignoring 
the large amount of shock heating, which increases c s , 
and considering adiabatic compression. Accordingly, for 
an n = 1, T = 2 fluid, P = KpQ, where k = 0.01 due 
to our pressure depletion. Fragmentation occurs when 
p ss p m 3, so Xj ~ 0.25 = 15AX. (Shock heating 
increases this coefficient.) Our resolution is then quite 
sufficient to resolve the Jeans length. 

We could not determine the final fate of this system. 
The four clumps may continue to collapse to black holes, 
or this collapse may be halted by heating-induced pres- 
sure. The system will certainly emit substantial amounts 
of gravity waves, both during the bounce and oscillation 
of the initially axisymmetric torus and during its rotation 
following fragmentation. To see this, we measured the 
gauge invariant Moncrief variables (or Zerilli func- 
tions) at the outer part of the grid [88j . We also measure 
the amplitudes of the two gravitational wave polariza- 
tions h+ and h x on the a;-axis at the edge of our grid. 
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FIG. 14: The gravitational wave amplitude h+, at a distance 
d from the source, for the 3D collapse, bounce, and fragmen- 
tation of star C. We compute h+ at the point (11.6 M, 0, 
0). 

Since the outer part of the grid is not in the wave zone, 
our measurements are only approximate. We find that 
the dominant mode of the emission is I = 2, m = 0, the 
quadrupole radiation generated by the axisymmetric col- 
lapse and bounce of the torus. The second largest modes, 
which are an order of magnitude smaller than the dom- 
inant mode, are I = 4, m = (octopole radiation from 
the axisymmetric collapse) and I = 4, m = ±4 (octopole 
radiation generated by the rotation of the four clumps) . 
In Fig. 1141 we plot h + on the x-axis, which contains con- 
tributions from all modes. The observed amplitude of 
this radiation from a star at a distance d from the Earth 
would be 

fc-lO-^V-JJ-V 1 (34) 

\M Q J ViooMpcy v ' 

The final evolution of this very interesting system can 
only be undertaken using a finer grid, presumably by 
employing adaptive mesh refinement (AMR) , and an im- 
proved shock-handling scheme in our code. 

VI. DISCUSSION AND CONCLUSIONS 

We have constructed a code to study the collapse of 
astrophysical objects to black holes by evolving the full 
coupled Einstein- hydrodynamics system in both 2+1 (ax- 



isymmetry) and 3+1 dimensions. When a black hole ap- 
pears, it is treated by introducing an excision boundary 
well inside the horizon. Our code is stable and conver- 
gent for all of the test problems and applications pre- 
sented here. As a test application, we study the collapse 
of rapidly rotating stars. Our conclusions regarding their 
ultimate fate agree with those of Nakamura and of 
Stark and Piran |8lj — namely, that spinning stars de- 
prived of their pressure support will collapse directly to 
black holes only if they are subKerr. This is the same be- 
havior observed for spinning configurations of collision- 
less matter j54[- We also were able to study the final 
state of the subKerr collapses by using our excision al- 
gorithm to extend the evolution far beyond what could 
be achieved without it. We find that even for a rapidly 
rotating star with q = 0.9, all the rest mass falls immedi- 
ately into the hole, with no disk formation, in agreement 
with Shibata [28|. For the case of supraKerr collapse, we 
found that the collapsing star hits a centrifugal barrier 
and bounces, forming a torus which fragments due to a 
nonaxisymmetric instability into four pieces. With our 
current computational resources, we were unable to de- 
termine the final fate of the four clumps. Systems like 
this one are sufficiently interesting as gravitational wave 
sources, that they should be pursued by further investi- 
gation with finer resolution, including AMR. 

Considering the stability of our excision algorithm over 
such a variety of applications, we believe that it has 
great promise as a tool for relativistic astrophysics in- 
volving the simultaneous presence of hydrodynamic mat- 
ter and black holes. Our current post-excision algorithm 
exhibits a gradual spurious decrease in total angular mo- 
mentum when applied at moderate resolution. However, 
this problem is not present in all coordinate systems (e.g. 
Kerr-Schild) and is reduced as the resolution is increased. 
We are currently investigating a number of ways to im- 
prove our algorithm and to apply it to other 2D and 3D 
problems of astrophysical interest. 
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